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Abstract 



Big Bang Nucleosynthesis (BBN) is studied within the framework of a 
two-parameter family of tensor-scalar theories of gravitation, with nonlinear 
scalar-matter coupling function a{tp) = ao + ao((p — (fo) + Tj(3(ip — (fo) 2 . We run 
a BBN code modified by tensor-scalar gravity, and impose that the theoreti- 
cally predicted BBN yields of Deuterium, Helium and Lithium lie within some 
conservative observational ranges. It is found that large initial values of a(ip) 
(corresponding to cosmological expansion rates, for temperatures higher than 
1 MeV, much larger than standard) are compatible with observed BBN yields. 
However, the BBN-inferred upper bound on the cosmological baryon density 
is insignificantly modified by considering tensor-scalar gravity. Taking into 
account the effect of e + e~ annihilation together with the subsequent effect of 
the matter-dominated era (which both tend to decouple ip from matter), we 
find that the present value of the scalar coupling, i.e. the present level of de- 
viation from Einstein's theory, must be, for compatibility with BBN, smaller 
than a 2 < 1(T 6 - 5 /T 1 (tt matter h 2 /0A5)~ 3 / 2 when > 0.5 . 
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I. INTRODUCTION 



The idea that Einstein's tensor gravitational field might be accompanied by a massless 
scalar partner was first suggested in the twenties by Kaluza 0. Since then, such tensor- 
scalar gravity theories have been studied in detail by many authors 0-0]. String theory has 
recently revived the motivation for considering gravitational-strength scalar fields, such as 
the (model-independent) dilaton or the (Kaluza-Klein-type) moduli (see, e.g., ||). In the 
simplest versions of tensor-scalar theories of gravity (those respecting the equivalence prin- 
ciple) the coupling between matter and the scalar field <p is described by a single "coupling 
function" a(ip), such that all physical mass scales get multiplied by a factor A((p) = exp [a (</?)] 
when measured in "Einstein units" (see below for the definition of the Einstein conformal 
frame). For instance, the Jordan-Fierz-Brans-Dicke (JFBD) theory is the one-parameter the- 
ory defined by a linear coupling function, a(ip) = ao (p. The JFBD theory is not an appealing 
alternative to general relativity because its only (dimensionless) parameter ao needs to be 
fine-tuned to a small value, < 10~ 3 , to be consistent with existing experimental data. 
[See, e.g., |§ for a discussion of the constraints on tensor-scalar gravity brought by solar- 
system, and binary pulsar, data.] By contrast, it has been shown that more general theories 
with nonlinear coupling functions, containing no small parameters, can be naturally compat- 
ible with experimental data because the cosmological evolution drives the background value 
of ip toward a value that minimizes the coupling function a(tp), thereby reducing by a large 
factor all the present observable effects of ip |IU|] , JTTJ. Within such cosmological-attractor 
models, the only regime in which the scalar field can play a quantitatively important role is 
early cosmology. 

Big Bang nucleosynthesis (BBN) is the earliest cosmological process that is physically well 
established and about which one has reasonably accurate observational data. Its importance 
for constraining many physical or astrophysical scenarios was first pointed out by Shvartsman 



12fl . For recent treatments see, e.g., ||13|- |i5|l . In particular, BBN is crucially used for 
deriving an upper bound on the cosmological baryon density 0& = Pbaryon/Pciosurc which, 
even for extreme parameters, is claimed to be Qb < 0.03 hr 2 < 0.2 |T3|-|T5|| (here, h = 
tfoAlOOkms -1 Mpc^ 1 )). One of the main motivations of the present paper is to examine 
whether the presence of a gravitational-strength scalar field ip, having a generic non-JFBD, 
i.e., nonlinear, coupling function a(<p), can significantly modify the standard limit on f^, 
while being still compatible with present gravitational experiments. A secondary motivation 
is to study the level of present scalar admixture to Einstein's gravity (observable, in principle, 
in precision tests of relativistic gravity) which is naturally compatible with BBN data, i.e. 
with observed light element abundances. 

Our analysis is more general and/or more exact than previous attacks on this problem. 
Indeed, we consider the case of a nonlinear coupling function a(ip) admitting a local mini- 
mum, while most previous studies restricted themselves to JFBD-type theories |TB|-|rg| (see 



also the related studies of BBN limits on the variability of the gravitational constant G |L2| , 
P"D|-P2"||). A recent work by Santiago et al. p|] has considered nonlinear coupling functions 
a(ip) of the type we shall study but their analysis is, in our opinion, less satisfactory than 
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ours in twoQ ways: (i) they approximated the effect of ip as a small, constant speed up factor 
£ n affecting only the value of the neutron/proton ratio at some effective "freeze-out", sup- 
posed to take place before and separately from e + e~ annihilation, without recomputing in 
tensor-scalar gravity the production of all the light elements, and (ii) they related the value 
of ip (and the corresponding £ n ) at freeze-out to its present value (po by integrating backward 
in time. Concerning the point i), we shall, instead, find that the beginning of e + e~ annihi- 
lation generates a source term for the evolution of (p which generically (especially for largish 
values of (3) cannot be neglected during freeze-out so that the speed up factor can vary a lot 
during the critical freezing of the n/p ratio. ( see Section 4 for a more detailed discussion). 
Concerning ii), the procedure of backward time-integration may generically lead to incorrect 
results because the f> evolution equation is similar to that of a damped harmonic oscillator 
ITJ . Integrating backward in time such an equation is equivalent to integrating forward 



in time an oscillator with negative friction, which introduces spurious runaway solutions. 
The origin of the introduction of these runaway solutions is simply that, when integrating 
backward in time, one lacks the correct "initial" condition for the present time derivative 
of <p, (po, corresponding to some <pa. [Indeed, <po is nonlocally defined by the full forward 
(damped) evolution of <p starting from the physically natural initial condition, (fi- m = 0, deep 
into the radiation era (see below).] We suppose that this pollution by runaway solutions is 
the root of the (incorrect) obtention of infinite peaks in the figures of Ref. [p3|| , corresponding 
to speed-up factors equal to one (which is physically forbidden, as we shall discuss below). 
[We also suppose that this pollution by runaway solutions, made more critical because of 
a direct numerical integration from the present time back to BBN time, invalidates the 
physical relevance of the work of Ref. ||24j| , and explains why their limits on the coupling 



parameter «q are about fourteen orders of magnitude smaller than the limits we get here.] 
However, we have not checked our supposition (concerning the effect of runaway solutions) 
by running reverse integrations ourselves. 

In this work, we shall avoid this potential problem of unphysical runaway solutions by 
integrating only forward in time. Another difference with previous work is that we shall 
compute, using a full BBN code modified by tensor-scalar gravity effects, the abundances 
of all the light elements: Deuterium, Helium 3, Helium 4 and Lithium 7. Indeed, though 
they are still large uncertainties on the primordial abundances of these elements, it is impor- 
tant to combine the predicted abundances of all the light elements with the corresponding 
observational bounds to get consistent constraints on tensor-scalar gravity. 

The present paper works within the framework defined by the following assumptions: (i) 
we consider tensor-scalar gravitation theories containing a single massless scalar field ip, (ii) 
the coupling function of ip is restricted to the two-parameter family a(<p) = a + ao((p — ipo) + 
\(3{(p — ipo) 2 (see Section II for the choice of this quadratic form), (iii) we require that deep 
into the radiation era the time derivative of (p vanishes (see our detailed discussion below), 
and (iv) we consider the case of a spatially flat Friedmann universe (the effect of non zero 
curvature has been investigated in Ref. |T0| ). 



x In addition, the solution for <p during the end of the radiation-dominated epoch is incorrectly 
given in Sec. IV D of Ref. [ 23 1 : e.g., their Eq. (4.37) should contain J\(x)/x and Y\(x)/x. 
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II. TENSOR-SCALAR GRAVITY THEORIES 



We consider tensor-scalar gravitation theories containing a single (massless) scalar field, 
assumed to couple to the trace of the energy-momentum tensor. This coupling means that 
the scalar source is (for bodies having negligible self gravity) proportional to the (inertial) 
mass, so that the equivalence principle is respected. The most general theory describing 
such a mass-coupled long-range scalar field contains one arbitrary coupling function J| 

A(<p) = exp(a(¥>)) ■ (2-1) 
The action defining the theory reads 

S = — ^- / d'xgl' 2 ^ - 2g^d^d v y) + S m ^ m , A a (<p)g^) . (2.2) 

Here, G> denotes a bare gravitational coupling constant, R* = g<£ v R*^ u the curvature scalar 
of the "Einstein metric" g*^ v describing the pure spin-2 excitations, and <p the long-range 
scalar field describing spin-0 excitations. [We use the signature — h ++ and the notation 
<7* = — det g* .] The last term in Eq. (|2.2j) denotes the action of matter, which is a functional 
of some matter variables (collectively denoted by ip m ) and of the (Jordan-Fierz) "physical 
metric" 

g„ v = A 2 (cp) 9 ; v . (2.3) 

Laboratory clocks and rods measure the metric g^ which, in the model considered here, is 
universally coupled to matter. In the BBN context, the standard laws of nongravitational 
physics (such as nuclear reaction rates and thermodynamical laws) will hold in their usual 
form when expressed in "physical units" , i.e. in units of the proper interval ds 2 = g^ v dx^ 1 dx v . 
For instance, Ncether's theorem applied to the matter action S m [i^ m )9tJ,u] yields the usual 
law of conservation of energy and momentum: 

V„T^ = 0, (2.4) 

where T^ v = 2g~ 1 / 2 5S m / 5g^ u is the material stress-energy tensor in physical units. In 
Eq. (|2.4|) the covariant derivative V„ is that defined by g^ v . 

On the other hand, the gravitational field equations of the theory are most simply for- 
mulated in terms of the pure-spin variables (g* , (p). From ( |2.2j ), they read 

r; u = 2 aw o vV + stt g* (t; v - 1 t* g ;)j , (2.5a) 

n g *(p = -47rG*a(^)T*, (2.5b) 

with T£ u = 2g* 5S m /5g* u denoting the material stress-energy tensor in Einstein units. 
It is related to the physical-units stress-energy tensor through 

g* va IT = = A\<p) T£ = A\ V ) g va f ^ . (2.6) 

The quantity a(tp) on the R.H.S. of Eq. ([2.5b]) plays a crucial role in the theory. It is the 
logarithmic derivative of the coupling function, 
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_ d In A(ip) _ dajip) 

<x(<P) = o = , (2.7) 

dip dip 

and it measures the basic (field-dependent) coupling strength between the scalar field and 
matter. The JFBD theory is defined by a linear field dependence of a(ip) = In A(ip) = ao ip, 
i.e. by a constant (field-independent) coupling strength a(ip) = ao [aft = (2u + 3) _1 in the 
notation of |4]]]. Generically, one might expect a nonlinear field dependence of a(ip) leading 



to a field-dependent coupling strength a(ip). It has been shown in Refs. [0], |25[ that all weak 
field ("post-Newtonian") deviations from general relativity (of any post-Newtonian order) 
can be expressed in terms of the values of a(ip) and of its successive (^-derivatives, starting 
with 

m^, (2. 8 ) 

at the present "vacuum expectation value" ipo of the field ip. Here ipo denotes the asymptotic 
value of ip at spatial infinity, at the present epoch. At the first post-Newtonian approxi- 
mation, deviations from general relativity are proportional to the well-known Eddington- 
Nordtvedt-Will parameters 

7 = Ttedd-l = -2ag/(l + ag), (2.9a) 
P = /3 Edd - 1 = +~ fa al/ {I + alf , (2.9b) 

where «o = a(y?o) and (3q = P(ipo). We see explicitly from Eqs. ( |2.9| ) that post-Newtonian 
deviations from general relativity tend to zero with ao at least as fast as a§. This holds true 
for weak-field deviations of arbitrary post-Newtonian order |25|| . By contrast, it was found 



in Ref. that strong-field deviations from general relativity do not tend to zero with ao if 



the parameter f3 = P(ipo) = d a(ip ) / d ip is sufficiently negative. However, if one considers 
(as we shall do here) the case where the cosmological-attractor mechanism of Refs. [ID 



jTT| takes place, the parameter (3 is necessarily positive. Indeed, Refs. QTOJ] , |TT| found that 
the spatial average of ip is attracted, during the cosmological evolution, toward a minimum 
of the coupling function a(ip). Therefore, the present cosmological value of ip, ipo = ip(to), 
is generically expected to be very close to a value ip m such that a m = d a(ip m ) / d ip m = 
and p m = d 2 a(ip m ) / d ip^ > 0. In such a case, all present deviations from general relativity 
(weak-field ones and strong-field ones alike) are expected to be very small, because they all 
contain a factor «q ~ [P m (ipo — ip m )} 2 "C 1- 

As was discussed in Refs. JTJJ|, [[FT], and is recalled in Section 3 below, the Einstein-time 



evolution for ip deep into the radiation dominated era can be approximated by ip+3H* ip ~ 0, 
which shows that ip decreases, during the expansion, as the inverse cube of the Einstein-frame 
scale factor (in usual cosmological parlance, one can say that ip contains only a decreasing 
mode). Therefore the natural initial conditions for ip deep into the radiation era are that up 
vanishes, while ip = ip- m has an arbitrary value. Starting from these values, ip will run down 
the coupling function a(ip) and will be attracted to the nearest minimum of a(ip). [Note that 
this shows, as emphasized in Ref. JTOJ, that the value of a(ip) irreversibly decreases during 



the expansion, so that the speed up factor must always be larger than one.] As the generic 
behaviour of a function near a minimum is parabolic, a(ip) = a(ip m ) + \ p m {ip — ip m ) 2 + 
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0((<f — fm) 3 ), it is plausible to assume that a{ip) is roughly parabolic in the whole interval 
around ip m containing ip in . We shall therefore follow Ref. || and study the paradigmatic case 
where the coupling function is quadratic in if, a quad (ip) = a Q + a Q (ip — ip ) + ~ /3 (ip — <fo) 2 . 
This case is the simplest generalization of the JFBD theory (gijfbd^) = a o + a o(v 9 — fo)) 
admitting a nonlinear coupling function, and susceptible of being cosmologically attracted 
to a point where p decouples from matter. By appropriately choosing the unit of length 
one can fix ao to zero. One can also conventionally fix the minimum of a((p) at (p m = 0. 
Therefore, without loss of generality, we shall write the general "quadratic model" in the 
simple form 

aq U ad(<^) = In v4 quad (v9) = - (3 ip 2 , (2.10) 

with (3 > 0. This model contains two free parameters: the quantity (3 = d 2 a((p)/d(p 2 
(= Po = Pm), and, either the initial value tp- iri of <p deep in the radiation dominated era, or, 
alternatively, the present value tpo of (p. The aim of the present paper is to investigate which 
values of (3 and ip- in (and consequently, (po) are compatible with BBN data. 



III. BIG BANG NUCLEOSYNTHESIS IN TENSOR-SCALAR GRAVITY 



As mentioned in Section 2, the standard laws of nongravitational physics will keep their 
usual form if one measures lengths and times in "physical units" , i.e. in units of the interval 
ds 2 = g^ v dx^ dx v = A 2 (ip) ds 2 . In a cosmological context, the physical interval will read 

ds 2 = -dt 2 + R 2 (t)d£ 2 , (3.1) 

where t is the physical cosmological time, and R the physical scale factor. Here 

dv^ 

di 2 = — + r 2 [dQ 2 + sin 2 9 d<j) 2 ) , (3.2) 

1 — kr z 

is the metric of a 3-space of constant curvature k — +1, or — 1. The corresponding physical 
expansion rate (Hubble parameter) is 

H = i,\nR. (3.3) 
dt 

The Einstein-frame counterparts of the above quantities are 

dsl = -dtl + Rl{U)dt 2 , (3.4) 



H* = ^-\nR*. (3.5) 



The relation ds 2 = A 2 (if) ds 2 gives the links 



dt = A(if) dt* , (3.6a) 
R — A(cp) R* . (3.6b) 
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Let us introduce the notation \ for the Einstein-time derivative of ip: 



at* dt 



(3.7) 



In terms of this notation, the links ( 3.6| ) give the following relation between the Hubble 
parameters, 



H=(A(<p))- l [H. + a(<p) x] 



(3.8) 



where a(<p) = da(<p)/d(p. 

We computed BBN by using Kawano's update p7] of Wagoner's code p8| . This code 
evolves thermodynamical variables and nuclear abundances as functions of the physical 
temperature T. All these physical evolution equations are unchanged in tensor-scalar gravity. 
The only modifications that need to be brought to this code are: (i) a modified expression 
for the physical expansion rate if as a function of the energy density, and (ii) two new 
first-order equations giving the evolution of the scalar field. Note that we assume three light 
neutrinos throughout this work. 

To get the modified expression for the expansion rate we need to write explicitly the 
modified Einstein equations ( 2.5a|) . In the present cosmological context they yield 

3 gP i?* 



R* dtl 



3 Hi 



k 

m 



47rG„(p* + 3p*) +2% 



8ir G* p* + x 2 ■ 



(3.9a) 
(3.9b) 



Here p* and are the Einstein-units energy density and pressure, linked to their physical 
counterparts by Eq. (|2.6|), i.e. by 



p* = A 4 p , = A 4 p . 



(3.10) 



where A = A(<p). As usual, the curvature term k/R^ in Eq. (|3.9b|) is negligible during BBN 
so that Eq. ( |3.9b|) gives the following result for if* in terms of the physical energy density 
and x = dip/dt* 



A*p + \x 2 



(3.11) 



The corresponding physical Hubble parameter is then computed by using Eq. 

It remains to write explicit evolution equations for the scalar field (p. Eq. ( |2.5b| ) yields 



I 3 H, ^ 
dtl * dt* 



— 47T G* a((p) (p* — 3p* 



(3.12) 



In terms of the variable x? Eq. ( |3.7| ), and of the physical time t, this gives the first-order 
system 
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^ = A-' X , (3.13a) 
at 

% = -A- l [3H, X + ^G,a(if)A A a], (3.13b) 
at 

where a denotes the following "source term" 

a = p-3p. (3.14) 

Around the period of primordial nucleosynthesis, the scalar source term can be approximated 
by 

5bBN - <Te = Pe ~ 3fi e , (3.15) 

where the index "e" denotes the contribution from electrons and positrons. Indeed, among 
the various contributions to the total energy entering Eq. (|3.11| ), 



P = P~, + Pe + Pv + Pm, (3.16) 

the massless (or, at least, ultrarelativistic) photon and neutrino contributions p 7 , p u satisfy 
o 1 = = 5v, while the non relativistic matter contribution p m (cold dark matter plus 
baryons) is negligibly small (it becomes, however, important later, around matter-radiation 
equivalence and during the matter dominated era). By contrast, it happens by coincidence 
that one of the most important phenomena of BBN, the freezing out of the neutron to 
proton ratio, takes place around a physical temperature Tp ~ 1 MeV, which is precisely 
when the e + e~ plasma starts becoming non relativistic before annihilating. In other words, 
the freeze-out of (n/p) takes place when there is a numerically large contribution a e to the 
source driving the evolution of the scalar field tp. This fact has been taken into account 
only in the coarse approximation of considering that the freeze-out takes place before, and 
separately from, e + e~ annihilation in Ref . p3f . Our full numerical integration of a scalar- 



tensor BBN code shows that a e drives a strong (generically oscillatory) evolution of (p at the 
same time (if not before) the (n/p) ratio freezes out (See Section IV C below). 

Summarizing, we completed an existing standard BBN code [E8J, |27J by adding: (i) 
Eqs. ( |3.8| ), ( p.llj ) for computing H in terms of p, <p and Xi an d (ii) Eqs. ( |3.13| ) for computing 
the t-evolution of tp and x- The source terms in these additional equations were taken to be 
( |3.15| ) and ( |3.16| ). We used for the crucial scalar source a e the expansion 



^ = ^^£(-r 1 ^^- (3.i7) 



7r2 _ "r : 



Here, K\ is a modified Bessel function, and z = m e /T (in units where h = c = k = 1). As 
this is an alternating series, it is numerically important to retain an even (and large enough) 
number of terms. We kept sixteen terms in the expansion fl3.17|) . [ By contrast, Kawano's 



code uses five terms. Note also that the coefficient in front of this expansion is misprintedF] 
in Ref. f23|.] See Eq. ( |3.23|) below for the exact, unexpanded expression of a e . 



We were informed by the authors of Ref. [23|, that they used the correct value in their code. 
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The rest of the code computes the evolution of thermodynamical variables and of nuclear 
abundances in terms of the physical temperature T. The link between physical temperature 
T and physical time t is given by the standard thermodynamical relation 

dT p + p ~ . 

-^ = -3-^— 3.18 
dt dp/dT 

The input data for each run of the code are: 

77 = (4/11) (n fe /n 7 ) in , y? in , Xin and (3. (3.19) 

The coefficient 4/11 is introduced so that 77 measures the final baryon-to- photon ratio, 
V = (nb/n 7 ) ou t, obtained after e + e~ annihilation. As is discussed below, when (3 is larger 
than about 0.2, the bare Newton constant entering Eqs. ( p. 11 ) and ( 3. 13b| ) can be 



taken to be numerically equal to the value presently measured in Cavendish experiments, 
G* ~ 6.672 x 10 -8 cm 3 g -1 s~ 2 . [This is because the cosmological attraction of ip toward 
(p m = is quite efficient when (3 > 0.2, so that the presently observable value G((fo) = 
G* A 2 (tpo) [1 + a 2 (ifo)] differs negligibly from G*.] The output data for each run are 



^out 



D/H , 3 He/H , F( 4 He) and 7 Li/H. (3.20) 



The ratio D/H = n(D)/n(H) denotes, for instance, the number of Deuterium nuclei per 
proton after BBN, while F( 4 He) ~ 4n( 4 He)/(ra(H) +4ra( 4 He)) is the abundance by weight of 
Helium 4. The other physical parameters needed in the code are taken to have their currently 
measured values (in physical units). In particular, the neutron (exponential) lifetime is taken 
to be r n = 888.54 s. 

An important issue (which, in our opinion, was not properly dealt with in previous work) 
concerns the initial values for ip and its derivative x = dtp/dt*. Our numerical calculations 
integrate the evolution for if in the normal direction, i.e. forward in time. If we insert 
our paradigmatic quadratic coupling function ( |2.10|) , leading to the simple linear coupling 
strength 

ttquadM = /3 if, (3.21) 

into Eq. ( |3.12j) , we see that (p(t*) satisfies the equation (here cp = d<p/dt*, cr* = p* — 3p*) 

<p + 3H*(p + 4:irG*a*(3(p = 0. (3.22) 

This is the equation of a damped harmonic oscillator with time-varying positive friction 
f = 3H*, and time- varying positive spring constant ^q(£*) = 47r (3 G*o*. Indeed, the scalar 
source cr* = p* — 3p* = A 4 (p — 3p) is always positive for a thermal distribution of particles 



g A m\ f 00 dqq 2 



p A - 3p A = / _ Ea = Jm\ + q 2 . (3.23) 

HA yA 2tt 2 Jo E A [exp(E A /T) ± 1] V A H K ' 



3 We do not mention 7 Be and 8 Li which we shall not use. 
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Here, tua is the mass of the considered particle, g& is its number of degrees of freedom, and 



the upper (lower) sign holds for Fermions (Bosons). As emphasized in Refs. [[nj, \\ 1] the 
positive damping in Eq. ( |3.22| ) means that when one starts the p evolution very early in the 
radiation-dominated era (where the total number of effective relativistic degrees of freedom 
is large) and/or for temperatures well away from any mass threshold T 7^ Sa (so that the 
'spring constant" term is relatively negligible in Eq. (|3.22 )) any initial "velocity" p- m = \ 



(satisfying Xm < V3H* in from Eq. ( |3.9b[ ) with k = 0) will be generically quickly (in a few e- 
folds) damped away (as said above, in these circumstances, up decreases like -R7 3 ). Therefore, 
a physically most natural requirement is to impose (as we do in this work) that the initial 
velocity Xm vanish when starting a BBN run at a temperature m M ^> T ^> m e . [The physical 
effects due to p would be stable under using only an inequality Xm *C H*, instead of a strict 
vanishing of Xm] Evidently, there exist mathematical solutions to Eq. ( |3.22j ) where <p/H* is 
not small before BBN. In view of the argument presented above, and in absence of a generic 
mechanism giving rise to large (i.e. very near a/3) initial values of <fi/H*, we think that such 
solutions are physically irrelevant. 

On the other hand, the initial "position" p- in is physically unrestricted, at least if the 
curvature parameter (3 is not very much larger than unity. [It was shown in Refs. ||10|| , [11 



that when (3 > 10, the effect of the previous mass thresholds (above the electron one) is very 
efficient in attracting <p toward zero, so that in this case one would expect tp- m before BBN 
to be small.] In this work, we shall formally leave pi n unrestricted, even when f3 ~ 100, to 
see what constraints BBN brings on this possibility. As the explicit numerical calculations 
we made start at the temperature T = 10 11 K ~ 10 MeV which is only 20 times larger 
than m e , we have slightly refined the initial conditions by taking into account the attracting 
effect of the source term a* in Eq. (|3.22|) , integrated over the temperatures > T ^> To- 



Using the asymptotic behaviour a e oc T~ 2 asf > m e , one finds that, starting from <p = <p in , 
X = Xm — at ray > T 3> T 3> m e , one ends up with the following corrected initial 
conditions at T = T : 

l p(f Q )^p ill --a m (^\ , (3.24a) 



4 \P, 




(3.24b) 

where ((pin). This "edge correction" refinement had, anyway, a negligible effect on our 

results. The important point to remember of this discussion is that the forward integration 
in time leads to a physically straightforward, and well motivated, choice for the initial time- 
derivative of (p. By contrast, any integration backward in time is equivalent to integrating 
an oscillator with negative friction (time-reverse of ( |3.22| ), with / = — / = — 3if*). Such 
an integration is unstable and may generically lead to including (without knowing it) a 
"runaway" mode in ip, i.e. a spurious backward-growing mode (equivalent to setting a large, 
non zero value of Xm at T 3> m e ). We think that this incorrect inclusion of spurious runaway 



modes affects the methods used in Refs. j23] and [24|, and may have invalidated some of 
their results. 
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IV. RESULTS OF NUMERICAL INTEGRATIONS AND THEIR COMPARISON 

WITH OBSERVATIONAL DATA 



A. Primordial abundances: observations 

The uncertainties in the "observed" primordial abundances of the light elements (Deu- 
terium, Helium 3, Helium 4, Lithium 7, . . .) are large and frought with systematic errors 
which are difficult to assess. See Refs. [[14]] for recent overviews of the observational 
situation. As the aim of the present paper is to set conservative limits on possible devia- 
tions from the "standard" (Einstein-gravity) BBN we wish to work with the most extreme 
(hopefully "3cr"-type) ranges for the possible "observed" primordial abundances. We shall 
take 



0.21 < F( 4 He) < 0.255, 
1.5 x 10~ 5 < D/H < 2.4 x 10' 
0.7 x 10~ 10 < 7 Li/H < 1 x 10 
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(4.1a) 
(4.1b) 
(4.1c) 



Let us comment briefly on our choices. The minimum value Y min = 0.21 of Eq. (|4.1a|) is 
the one admitted (as extreme possibility) by Refs. |T3[ and |14]], while Y max = 0.255 is an 
extreme maximum value suggested by Ref. P9[| . The minimum value for D/H, Eq. Q4.1b|) , 
is from [|I4J], while the maximum one was derived by assuming that the "high" Deuterium 
value (with its one-sigma error bar) reported by Ref. j30| might be correct. We are aware 
of the observational results, and arguments, of Tytler and collaborators converging toward 



a "low" Deuterium value, say |]31 | 



D/H = (3.4 ±0.3) x 10' 



(4.2) 



but we wish to be very conservative to draw secure limits on tensor-scalar gravity. We shall, 
however, comment below on the effect that a more precise knowledge of the primordial 
Deuterium abundance, of the type of Eq. ( |4.2|) , would have on our results. The minimum 
value for Lithium 7 is taken from Ref. |T4|], while the maximum one is the extreme value 
chosen in Ref. 



la 



Finally, we have also computed the BBN yield of Helium 3. We used 
an extreme range on the (observationally better controlled) sum (D + 3 He) /H. However, 
we found that the corresponding limits do not give constraints going beyond those deduced 
from the BBN yields of Helium 4, Deuterium and Lithium. 



B. Small f3 case (constant speed up factor) 

Let us first consider the case where the curvature parameter j3 of the quadratic coupling 
function (|2.10| ) is small, say f3 < 0.2. In the limit f3 — > 0, a qua d(y2) = (3ip — > 0, and the 
evolution equation ( p\22j ) for ip becomes simply <p + 3f/* = 0. The generic solution ip 
of this equation tends to an arbitrary constant ip = C as time evolves in the (normal) 
positive direction. [The other solution is the decreasing mode <pd ec = D/B% which is quickly 



damped away] It was shown in Ref. [10] that such a behaviour, ip ~ constant, is a good 



approximation to what happens during the e + e annihilation era when (3 is small. Indeed, 
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Eqs. (4.20), (4.21) of Ref. |T(J show that the fractional change of <p across the e + e threshold 



is of order A e ip/ip ~ —k e f3 ~ —0.16/3, where the smallish coefficient k e is proportional to 
the ratio of the number of helicity states of an electron to the total number of relativistic 
degrees of freedom when T > m e . If we limit ourselves to values of (3 < 0.2, the fractional 
change of ip will be less than 3% and can be neglected. This means that, when (3 < 0.2, the 
effect of tensor-scalar gravity on BBN is simply obtained by replacing ip by a constant ip- m 
and x — <P by in Eqs. ( ft.8| ) and ( |3.11| ) giving the physical expansion rate. This yields 



^(small/3) = y ^^A 2 (ip in ) p. (4.3) 

By comparison, the standard value of the expansion rate in general relativity would read 

^standard = ^fe Q N p, (4.4) 

where Gm is the presently observed value of Newton's constant (as measured in Cavendish- 
type experiments). The value of this observable quantity in tensor-scalar gravity is 

G N = G = G,A 2 (p, )[l + a 2 ]. (4.5) 

Therefore, the small-/? limit of the tensor-scalar expansion rate is obtained from the standard 



value by the following rescaling of the standard gravitational constant [10 



= GJ 2 ( fa ) A 2 (<p in ) 
9 ~ Go A*{<po)[l + afl- { ' } 



Equivalently, one could say that the expansion rate ( |4.3|) differs from its standard counterpart 
by a constant "speed up" factor £ = H/H st3Ildaid = ^fg. In our quadratic model (|2.10| ), the 
natural logarithm of g reads, in terms of a in = a(<p- m ) = {3 <pt 



In g = - (a 2 m - ajj) - ln(l + ag) ~ - (a 2 n - a 2 ) - a 2 . (4.7) 

Here, we used the fact that Oq < 1, as known observationally, and also consistently deduced 
below from BBN data. To express lng in terms of ao and (3 we need to relate ao and a- m . In 
our present approximation, a- in = (3 <p- in ~ a ou t = (3 (p on t, where (p out is the value of ip at the 
end of nucleosynthesis. The link between (p ont and the present value ipo of the scalar field 
is obtained by integrating Eq. ( |3.22j ) (with initial conditions <p = (font, = 0) from the end 
of nucleosynthesis up till the present time. The evolution of <p during this time span comes 
from its coupling to the total matter (baryons plus cold dark matter). The corresponding 
scalar source a* m = p* m — 3p* m ~ p* m has a small effect on the evolution of <p during the 
radiation era but becomes very important during the subsequent matter-dominated era. It 
was very generally shown in Ref. |H| that the £*-time evolution equation (|3.12|) for ip (which 
must a priori be solved together with the evolution equation for the scale factor R*) can be 
conveniently rewritten as a decoupled equation for the evolution of <p with respect to the 
parameter 
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V 



H* dt* = In i?* + est 



(4.8) 



which measures the number of e-folds in the Einstein frame. The latter p-time equation was 
approximately solved in Ref. [TI| for the present case of the transition between the radiation 
era and the matter era. For simplicity, and because it is the theoretically most plausible 
value, we assume here and below that the spatial curvature k = 0. [ The effect of nonzero 
curvature has been investigated in Ref. ||10|| .1 The result for small (3 reads 



Q'o 



y^out ciout 



1 + r 
"27 



exp 



4 



(1 - r)p 



-Ppo 



Here, r 



3/3/3 = 1 - 4/3/3 + 0(/3 2 ), and 

po = in(;w#: quivalence ) 



(4.9) 



(4.10) 



denotes the number of Einstein-frame e-folds separating the present moment from the time 
of equivalence between matter and radiation: 1 = (p* m /p* r ad) eqmvalence = (p m /p ra d) eqmvalence - 
The corresponding physical redshift is 



Rn 



70 



P 



^equivalence 



2.404 x 10 4 Q m h 2 



(4.11) 



Here, Q m = p^/p^ is the present ratio of the total matter density to the (physical-units) 
"closure density" ft = 3H$/(&rrG Q ) and h = H Q /(100 kms" 1 Mpc" 1 ) is the reduced Hubble 
parameter. Note that the product 



~o 

Pm 



1.879 x 10- 29 gcm- 



(4.12) 



is independent from the value of H and depends only on the actual value of the total matter 
density. Finally, using the link ( |3.6b| ) between R and i?*, the number of Einstein-frame e- 
folds (fliy) reads fH| 



Po = \nZ + a(ip. 



out/ 



a(<p 



which gives numerically 



po = 10.09 + ln(O m h 2 ) + a{ip c 



a(<Pc 



(4.13) 



8.19 + ln(O m h 2 /0.15) + a(ip out ) - a(ip ) . 



(4.14) 



Here, we introduced as fiducial value for the matter density the value (Q m fi 2 ) Muaal = 0.15 
which corresponds roughly to the currently favoured values Q m ~ 0.3, h ~ 0.7 (it could also 
correspond to Q m ~ 1, h ~ 0.4). 

Using the link ( |4.9|) we can express a in , in Eq. ( |4.7|) , in terms of a : «i n — «out : 
This gives 



h/3po 



\ nfJ ~ !^[(?Pp°-1-/3\ 



(4.15) 
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Note that, though we are considering here smallish values of (3 ~ 0.1, we are not entitled 
to expanding the exponential in Eq. ( [4.15|) . Indeed, we see from Eq. ( 4.14 ) that 2p ~ 16 



is rather large. When (3 is very small (with 2p (3 <C 1), the expansion of the R.H.S. of 
Eq. ( [4.1 5|) gives a finite limit 



(ln<7V^ = (2po-lK. (4.16) 

It is easily checked that the R.H.S. of Eq. ( |4.15| ) is positive (because p > 1/2). Therefore, 
we conclude that, in the limit of small f3, tensor-scalar gravity leads to speeding up the 
expansion rate: £ = ^fg > 1. It is, however, interesting to leave first unrestricted the sign 
of Ing and to study the effect on BBN of replacing Gn —> g Gn for both signs of In g. We 
have run BBN codes with varying values of g and rj = fib/fi 7 , and computed the resulting 
abundances of light elements, comparing them to the extreme observational ranges given in 
Eqs. ( |4.1|) . The results of this comparison are given in Fig. [I] which represents the allowed 
regions in the (rj, g) plane. The total phenomenologically allowed range for g would be 

gmm < g < sw with 

£ min ~ 0.64 , g max ~ 1.27. (4.17) 

In our case g is restricted to be larger than one and we have the firm constraint 1 < g < 
5Wx- Note that the maximum allowed value of g is constrained by the upper bound on the 
deuterium, combined with the upper bound on Helium 4. To illustrate the effect of choosing 
a much smaller range for deuterium, we have also plotted on Fig. [I] (in dashed lines) the 
limits on g obtained by using as confidence interval the (la) one suggested by the recent 
results of Tytler and collaborators ||31|| , Eq. ( |4.2j ). This would give the stronger constraint 

S'max — 1-1. 

Note also that the shapes of the curves bounding the maximally allowed region are such 
that, when g is only allowed to vary above one, the corresponding allowed interval for the 
values of the baryon to photon ratio i]i = 10 10 x rj does not change significantly: roughly 
1.5 < r] W < 8.5. Therefore, when (3 < 0.2, the addition of a scalar component to gravity 
does not change the upper bound on fi& h 2 deduced from BBN. The case (3 > 0.2 will be 
considered below. 

Let us now turn to the constraints on the parameters «o and (3 of tensor-scalar gravity 
following from BBN. Using the result Q4.15D , the constraint 1 < g < g max translates in the 
following bound on 



a l< e 2 Ppo t l _ f3 ln ^x- (4.18) 
In particular, the limit ( |4.16j ) when (3 <C l/(2p ) yields 

(«o)/3^0 < - hlffmax, (4.19) 

which is numerically of order («o)/3^o ~ 0.015. Note that, given a value of Q m h 2 , the precise 
value of po defined by Eq. ( |4.14| ) must be obtained by iteration because the small (but non 
negligible) additional term 



14 



depends both on ao and on p . However, in our present approximation of small (3, we can 
replace flfclgp in ( ^20l) to get 



/ \ maximum 

I, "out a 0jsmall/3 



e 2/3po _ 1 



2 ( e 2/3po -1-13) 



In <7max — 2 ^ S'max 



(4.21] 



We then deduce from this result that the value of po which can be used for estimating the 
maximum BBN-allowed a 2 , is 

(po)smafl/9 ^ 8.19 + ln(Q m /i 2 /0.15) + - In # max 



8.31 + ln(fi m /i 2 /0.15). 



(4.22) 



C. P > 0.2 case 

Let us now consider the complementary case where (3 > 0.2. In this case, the total 
attraction due to the matter era is quite large. Indeed, the ratio ®.l/a% ut decreases fast 
with P, and even for /? = 0.2, one finds from Eq. ( |4.9| ) ao/a 2 ut — e~ 2/3po ~ 0.04, where 
we used po ~ 8. Therefore, when /3 > 0.2 we can, as was mentioned in Section 3, neglect 
the difference between G* and Go = G* A 2 (<p )[l + a 2 ]. This allows one to run (in the 
forward time direction) the tensor-scalar- modified BBN code described in Section 3, with 
the numerical value ~ 6.672 x 10~ 8 cm 3 g _1 s~ 2 . Let us recall that the input parameters 
which must be chosen in each run are: /3, the baryon-to-photon ratio 7710 = 10 10 x r], and 
the value f>- m of the scalar field before freeze-out and e + e~ annihilation. Instead of working 
with (p in , it is physically more appropriate to work with the corresponding value of the 
(logarithmic) coupling function a in = a(y? m ) = f P <^ 2 n - For each value of the triplet (P, r], a in ), 
the code then computes the abundances of light elements produced by the BBN, and the 
value ip ont of the scalar field after the end of nucleosynthesis (together with the corresponding 
flout = a(y?out) and a out = cu((p ou t)). Let us recall that for temperatures well above and well 
below the electron-annihilation threshold T ~ m e , the source term a e on the R.H.S. of the 
evolution equation (|3.13b| ) for (p becomes negligible, while the friction term remains, so that 
if stops evolving. [This approximation is correct as long as T < m M and T 3> T equivalence .] 
On the other hand, when T is between 5m e ~ 2.5 MeV and 0.2 m e ~ 0.1 MeV, which is 
a crucial period during which weak interactions freeze (the (n/p) ratio freezes out around 
T ~ 0.8 MeV), and nuclear reactions start to build up light elements, the source term a e 
is numerically important and causes if to oscillate around (p — 0. As discussed in detail in 
Ref. [|lTJ| these oscillations are very vigorous when P ^> 1. Therefore, we have a physically 



non trivial situation where the physical cosmological expansion rate H (given in terms of 
ip by Eqs. ( |3.8|) and ( |3.11| )) varies in a complicated manner precisely during the crucial 



stages of BBN. This situation cannot be analytically approximated, and must be tackled 
numerically. We can visualize the results of our numerical simulations by representing, for 



15 



each value of /?, the contour levels of the light element abundances in the (n, a- m ) plane. The 
topology of these contour levels depend on how large (3 is. 

When (3 < 16, the contour levels look like continuously deformed versions of the small-/? 
case, i.e. the case of a constant speed up factor £ = yjg shown above in Fig. [l| For instance, 
the case /3 — 10 is shown in Fig. and should be compared with the part of Fig. [I] above 
the g — 1 line. Therefore, when (3 ^ 16 the upper bound on the value of a- m is given, 
as in the small-/? case, by the upper left corner of the "triangle" defined by the (D/H) max 
and the He max lines. We have checked that there is approximate continuity between the 
small-/? case and the (3 > 0.2 case: Indeed, by running a BBN code for (3 = 0.2 we find a 
maximum allowed value a™ ax = 0.125 (together with a corresponding a™ ax ~ 0.109) which 
approximately agrees with ~ In g max ~ 0.12, for the g mSuX ~ 1.27 obtained above in the small- 
(3 approximation. [Here, we used Eq. ( |4.7| ), in which we neglected the fractionally small «q 
terms.] As in the small-/? case, we find that the shapes of the curves bounding the allowed 
region are such that the corresponding allowed interval for 7710 does not change significantly. 

When /? > 16, the "triangle" of Fig. § disappears because many of the contour lines tend 
to straighten out and become parallel. This situation is illustrated in Fig. |3| for (3 = 40. 
This case illustrates two interesting phenomena. First, though the shapes of the contours 
have changed a lot, they are still such that the allowed interval for 7710 does not change 
significantly: roughly, 1.5 < 7710 < 8.5. Second, Fig. |] shows that the initial value ai Q , i.e. 
the initial speed-up factor £ in = A in / (A \Jl + a^) ~ A- m = exp(a in ) can be extremely large 
(for instance, the largest value a- in ~ 5 of Fig. |3] corresponds to £i n ~ 150, and much larger 
values are possible) and still lead to compatibility with the observed abundances of light 
elements. This result contradicts the usual expectation that the BBN data are tight enough 
to constrain to a small level any deviation from standard gravity around the epoch of BBN. 
Here, we have an example where, just before the BBN starts, the tensor-scalar-predicted 
expansion rate is several orders of magnitude larger than the Einstein-predicted one, and 
where the BBN test is still passed. This paradoxical result raises two questions: (i) How can 
an abnormally large expansion rate, around the (nominal) time of freeze-out of the neutron 
to proton ratio, lead to quasi- normal light element abundances ?; and, (ii) Does this imply 
that, when /? > 16, the present value of a% (which measures the deviation between general 
relativity and tensor-scalar gravity) can be large? 

To answer the first question, we plot in Fig. [| the coupling function a(ip), which is 
essentially the logarithm of the speed-up factor, as a function of the inverse (physical) 
temperature 1/T (we use as input values (3 = 40, and a in = 2). We plot on the same 
Figure the evolution with temperature of: the actual neutron to proton ratio in tensor- 
scalar theory (solid line), the (n/p) ratio in the standard general relativistic case (dashed 
line), and the instantaneous equilibrium value for this ratio: (n/p) eq = exp(— Q/T), where 
Q = m n — m p ~ 1.293 MeV (dotted line). If we define (as is often done) a nominal freeze-out 
temperature as the temperature T* where the actual (n/p) ratio starts to deviate significantly 
from (n/p) eq (a sign that weak interactions become slow compared to the expansion rate), 
we see that the value of T* in tensor-scalar gravity is larger than in general relativity, as 
expected from the fact that the expansion rate is much larger (a(ip) being initially frozen at 
a ~ a in = 2). If that were all, this would lead to a significantly larger value for the (n/p) 
ratio. However, Fig. |4] illustrates two physical phenomena that modify this naive conclusion. 
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First, around T ~ 0.8 MeV, the source term a e starts to have a strong effect on a(<p): it 
makes a(ip) drop precipitously toward zero, where it bounces twice to end up being frozen 
(when T < 1/20 MeV) at a very small limiting value a out — 0.016. The second important 
phenomenon illustrated by Fig. ^ is the fact that the freeze-out phenomenon is not something 
which takes place at a precise moment in time, but is an integrated phenomenon. It is clear 
from the Figure that the approximate plateau reached by (n/p) (before it drops around 
1/15 MeV when nuclear reactions start) is due to the integrated effect of an evolution which 
takes place before, during and after the precipitous fall of a(ip) toward small values. It is 
therefore impossible to estimate analytically the final "plateau" value of (n/p) from the sole 
knowledge of the nominal freeze-out temperature. In particular, we see that (n/p) is initially 
larger than the standard GR value, but ends up, in this example, being somewhat smaller. 
This Figure illustrates the necessity of resorting to a full numerical BBN code. 

Let us now address the second question: What does an initially large value of a- in = a(ipi n ) 
imply for the present value of o/q ? Can q.q be large ? The answer is no, for the following 
reason: when ft is large, the source a e , active around T ~ m e , is very efficient in attracting 
ip toward zero (through a damped-oscillatory evolution). Moreover this attraction effect is 
amplified by nonlinearities so that, for a fixed ft, the post-BBN value a ont = a(ip ont ) does 
not increase monotonically with a- iri = a(y>j n ), but rather decreases for large a in values after 
having reached a maximum when a- in ~ 1. This behaviour is illustrated in Fig. [5[ We see 
from this Figure that, in the plane of Fig. |^, a out (considered as a function of a; n ) reaches 
the maximum value a™^ x ~ 0.028 when a- m ~ 0.7. 

For some particular values of ft (near 23, 59, . . .) the curve giving a out versus a in develops a 
second maximum which becomes higher than the first (like in a first-order phase transition). 
This is illustrated in Fig. |6] which shows the evolution of the curves between ft = 22, 23 and 
24, and between ft = 58, 59 and 60. In both cases, the formerly unique maximum located 
near the origin gets superseded, roughly when its location along the a; n axis decreases down 
to ~ 0.2, by a secondary maximum (which develops abruptly for ft > 21, or ft > 56) located 
around a in ~ 1.7. As ft further increases, the first maximum gets very small and disappears, 
while the position, along the of the second maximum continuously slides toward 

the left (so that the "first maximum" in the bottom panel of Fig. 5 is the evolved form of 
the second maximum in the top panel). 

In all cases, the maximum possible value of a out , or better of the difference a out — ao 
(which matters most when ft is small), i.e. the maximum deviation from general relativity 
after the end of nucleosynthesis, compatible with BBN data, is quite small. This is shown 
in Fig. [7] which plots the maximum allowed a out — ao in function of ft. When ft > 0.2, the 
maximum value of a out — a ~ a out is obtained from the results of the present sub-section, 



while for ft < 0.2, a out — a is obtained from the small-/3 result, Eq. ([4.21 



V. BBN LIMITS ON THE PRESENT DEVIATIONS FROM EINSTEIN'S 

THEORY 

In the previous Section, we have shown what limit BBN data put on a out = a(ip ont ) as a 
function of ft. It remains to translate this limit on the value of tp at the end of nucleosynthesis 
into a limit on the present scalar-coupling parameter ao = a(tpo). To do this we need to 



17 



integrate the ^-evolution through the transition between the radiation era and the matter 
era, up till now. As we said above, ip satisfies a decoupled equation in terms of the parameter 
p of Eq. (}4.8|) |T0|. In the present case of the transition between the radiation era and the 



matter era, it was further shown in Ref. JTTJ] that the latter p-time decoupled evolution 



equation can be approximately rewritten as the following hypergeometric equation 

x(x + l)%<p+ + d x <p + ^/3<p = 0. (5.1) 

Here x = e p = R*(t*) / Rl qmvalence denotes the (Einstein frame) redshift separating the current 
time U from the time of equivalence (when p* m (t^ uivalence ) = P* radiation (t^ uivalence )). The 
initial conditions for tp, deep into the radiation era (i.e. for x — > 0), mentioned in Section 
IVB above, select uniquely the solution ip(x) = <p ou t F[a, b, c; —x\. Here, F[a, b, c; z] denotes 
the usual (Gauss) hypergeometric series. The values of the parameters are 

3 3 /3 / 3\~ 

a = --iu , b = - + iu 1 c = 2 ; = W - f /3 — -J . (5.2) 

Finally the ratio between the present value of <p and its value at the end of nucleosynthesis, 

(p a 



"Pout «out 



= F m (P,Z, Q ) (5.3) 



is given by 



F m = F[a, b, c; -Z*] ~ ^ 2 2iw I&^l e -f *> j*t* + <_> _ iu} ) . (5.4) 



0r r (2iu + 



In the last form (made of two terms obtained by changing the sign of u) we have used the 
asymptotic behaviour, for large argument, of the hypergeometric function. Here, Z*o = e Po 
denotes the Einstein-frame redshift separating the present moment from equivalence. Its 



natural logarithm p = In Z* is given by Eq. ( 4.14 ). The result ( |5.4| ) is valid both when (3 < 



3/8 (in which case uj is pure imaginary, with iu = y | (| — /?)) and when (3 > 3/8 (u real). 

It is easily checked that when (5 is small the result Q5.4Q agrees with the more approximate 
Eq. ( |4.9|) . Using Eq. ( |5.3| ) we can compute 0% in terms of a out = \ (3 <p 2 nt = a out/ (2/5) 

a 2 Q = F%(J3, Z #0 ) a 2 oat = 2/3F%(j3, Z* ) a ont . (5.5) 

Using Eq. (|5.5|) we can, for a given value of Z*o, translate the BBN limit of Fig. [7] on 
a out — O'O into a limit on the present scalar-coupling parameter Oq. As before, when f3 > 0.2, 
a out — — o out and the limit on «q is obtained from Eq. (|5.5|) , while, when (3 < 0.2, one 
uses the small-/? limit, Eq. ( |4.18j ). This limit is represented in Fig. [8] (solid line) for our 
fiducial value Q m h 2 = 0.15, in Eq. ( |4.14|) . The corresponding limit on the product (3 a 2 , is 
represented in Fig. |9|. We recall from Eqs. ( |2.9| ) that the quantities represented in Figs. || and 
[| are directly related to the parameters measuring post-Newtonian deviations from general 
relativity. As ap < 1, we have 7 = 7Edd — 1 — — 2 ckq, and (3 = /?Edd — 1 — 2 P a o- 
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For comparison let us recall the current direct experimental limits on a 2 , and 8 oiq- The 



measurement of the Shapiro time delay by the Viking mission [ 32|| , as well as some VLBI 



measurements [Q, yield the lex bound |7| < 2 x 10 3 which translates into «q < 10 3 . A 
more stringent limit follows from the Lunar Laser Ranging experiment [[34] which yields the 
la bound -1.7 x 10~ 3 < 4/3 - 7 < 3 x 10" 4 . This translates into 

- 8.5 x 10~ 4 < (1 + A,) al < 1.5 x 10~ 4 (la) . (5.6) 

In the framework of the present paper, the parameter /3q = 8 is restricted to be positive. 
Therefore, the most stringent limit is the right inequality in Eq. ( |5.6|) , i.e. 

9 1.5 x 10~ 4 „ 9 8 d , > 

< ~/3~+i ' K 1,5 x 1 ' ( 5 - 7 ) 

[Strong-field tests put also limits on a 2 , but they are slightly less stringent than the above 
when 8 is positive ||.] We have indicated on Figs. § and |9] the empirical bounds ( |5.7j ) in 
dashed lines. 

Figs. [8] and |9| show that BBN data put, as soon as 8 > 0.3, more stringent limits on the 
possible post-Newtonian deviations from general relativity than present experimental data. 
For our fiducial matter density VL m h 2 = 0.15 the level of deviation compatible with BBN 
data is (for 8 ^ 0.5) constrained to the level 8 < 10~ 6 ' 5 . As discussed in Ref. |10| the 



attraction factor F^(8, Z^) in Eq. QJT^) scales (for 8 > 3/8) like (fi m h 2 )~ 3 / 2 . Note that the 
corresponding deviation level a 2 , < 10 — 6 5 /3 — 1 (fl m /i 2 /0.15)~ 3//2 gets significantly smaller if 
fl m h 2 turns out to be of order one. Note finally that the attraction factor F m is dominated 
by the evolution during the matter-dominated era and is only negligibly modified (in view 
of our other approximations) if one takes into account the fact that our recent cosmological 



evolution may have been curvature-dominated |10[ or A-dominated. 



VI. CONCLUSIONS 

The most natural, and best motivated, alternatives to general relativity are the tensor- 
scalar theories of gravitation. We have considered the class of tensor-scalar theories which 
are dynamically attracted, during their cosmological evolution, toward general relativity. 
The paradigmatic example of this class is defined by the quadratic coupling function a(ip) = 
\8 (p 2 with positive 8- We studied Big Bang Nucleosynthesis in this quadratic model as a 
function of 8, of the initial value ip- in (at temperatures ~ 10 MeV) of the scalar field, and of 
the baryon-to-photon ratio rj. [We assume three light neutrinos throughout this work.] We 
imposed that the theoretically predicted BBN yields of light elements (Deuterium, Helium 4 
and Lithium 7) be compatible with some conservative ranges for the corresponding observed 
primordial abundances, Eqs. (f4TT|). 

Our first conclusion is that the BBN-inferred upper bound on the cosmological baryon 
density Q b h 2 = 3.66 x 10 7 (T /2.726 K) 3 7/ is quite robust under the addition of a scalar 
component to gravity. The standard bounds on r/10 = 10 10 rj, namely 1.5 < r/10 < 8.5, 
corresponding to 0.0055 < Q^h 2 < 0.031, are insignificantly modified even when considering 
large initial values of the scalar field, corresponding to an initial cosmological expansion 
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rate much larger than its standard general relativistic value. This is strikingly illustrated in 
Fig.! 

Our second conclusion is that, even in the cases where, before BBN, tensor-scalar gravity 
is very different from Einstein's gravity, the presently observable deviations from general 
relativity are constrained by BBN to be quite small. The BBN-limits on the two independent 
weak-field "post-Einstein" parameters a 2 , and f3 a 2 , (linked to the usual post-Newtonian 
parameters through 7Edd — 1 — — 2a 2 ,, /^Edd — 1 — \ (3 ck 2 ,) are exhibited in Fig. |8] and Fig. |S|, 
in the fiducial case where the total matter density is £l m h 2 = 0.15. The BBN-limits in 
Fig. | and g scale (for (3 > 3/8) like (Q m /i 2 /0.15)" 3/2 . The dashed lines in Figs. 1, 2 
and 3 show also that if, in the future, one gets more stringent bounds on the Deuterium 
primordial abundance this will tighten the BBN-limits on a 2 , for (3 ^ 16, without affecting 
them for (3 ^ 16. As soon as (3 ^ 0.3 the BBN-limits we obtain on possible deviations from 
Einstein's theory are much stronger than the present observational limits from solar-system 
or binary-pulsar gravitational experiments. They provide motivations for experiments which 
push beyond the present empirical bounds on the basic scalar coupling parameter a 2 ,. 
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FIGURES 




^10 



FIG. 1. Level contours of BBN yields as functions of the baryon-to-photon ratio 
r/io = 10 10 77 = 10 10 nf,/n 7 and of a (squared) speed-up factor g = £ 2 = Q eSectlve j 'G?Newton- The 
solid lines correspond to the conservative observational bounds of Eqs. ( |4.1| ). The allowed region 
is the inside of the curved parallelogram defined by the He and D/H lines. The dashed lines, given 
for illustration, correspond to the Deuterium la range of Eq. (|4.2|), 
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|8 = 10. 




FIG. 2. Level contours of BBN yields when (3 = 10 as functions of the baryon-to-photon ratio 
7710 and of the initial value of the scalar coupling function a- in = a(y?i n ). The lines are defined as in 
Fig. ||. The allowed region for the conservative bounds is the truncated quasi-triangle defined by 
He max , (D/H) max , (D/H) m j n and the horizontal axis. 
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|8 = 40. 
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FIG. 3. Same as Fig. || for (3 = 40. The allowed region has now opened into a vertical strip 
essentially defined by the D/H lines. The wiggles in the He m j n contour signal the presence of some 
numerical noise. 
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Coupling function and (n/p) ratios vs. temperature 
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FIG. 4. Evolution of the coupling function a(ip), and of the neutron to proton ratio, in function 
of the inverse of the physical temperature (for [3 = 40 and a m = 2). The solid line is (n/p) in 
tensor-scalar gravity, the dashed line is (n/p) in general relativity, while the dotted line is the 
equilibrium value of the ratio (n/p). a out indicates the asymptotic level reached by a(<p) for small 
temperatures. 
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FIG. 5. Value of the final scalar coupling function a ou t = o-ifout), after BBN, in function of 
its initial value a; n = a((f[ n ) for (3 = 40. Because of the presence of a maximum in this curve the 
unlimited allowed region of Fig|3| still leads to a tight bound on a out . The wiggles for a; n > 3 are 
due to numerical noise. 
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a out 




FIG. 6. Same as Fig. | for (3 = 22,23 and 24 (top panel), and for (3 = 58,59 and 60 (bottom 
panel). This figure illustrates the competition between two maxima in the curve a ou t(oin)- F° r 
/? > 23 and > 59 it is the second maximum which determines the maximum allowed a ou t- 
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FIG. 7. Maximum BBN-allowed value of the difference a out — ao as a function of (3. 
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FIG. 8. The solid line represents the maximum BBN-allowed value of the present 
scalar-coupling parameter function of j3. For comparison the dashed line represents the 

most stringent direct empirical bound on Oq. 
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